Boundary Layers and Singular Perturbation 

Lectures 16 and 17 

Boundary Layers and Singular Perturbation 

A Regular Perturbation 

In some physical problems, the solution is dependent on a parameter e. When the parameter e is 
very small, it is natural to expect that the solution not be very different from the one with e set to 
zero. Consider for example the motion of a harmonic oscillator of unit mass with the spring constant 
equal to (1 + ee~'). The equation of motion for this problem is 

x + (1 +ee~')x = 0, (9.1) 

with 

*(0,e) = 1, i(0,e) = 0, (9.2) 

where x(t, e) is the displacement of the harmonic oscillator and x is the time derivative of x. 

Equation (9.1) can be solved in a closed form and the exact solution is a Bessel function. But if e 
is very small, the spring constant is approximately equal to unity at all times t > 0, and to the 
lowest-order of approximation, we expect that we may ignore its small deviation from unity. Thus 
we have 

x(t,e) « x(t,0) = cost. (9.3) 

To find the small correction due to the small deviation of the spring constant from unity, one may 
anticipate that, since the deviation is of the order of e, this small correction is of the order of e as 
well. Indeed, since all terms in (9.1) are analytic functions of e, we will try to see if there is a 
solution in the form of a Taylor series of e: 

x(t,e) « x (t) + exi(t) + ■■■ + e n x„(f) + ■■■, (9.4) 

where x„(t) is independent of e. Substituting (9.4) into (9.1) and (9.2), we get 

+ l)[*o(f) + e*i(f) + ■■■] + ee-'[x (t) + €x l (t) + •••]= 0, (9.5) 
*o(0) + eti(0) + - = 1, (9.6) 

and 

*o(0) +e*i(0) + ••• = 0. (9.7) 
Setting 6 in the equations above to zero, we get 

x'o + xo = 0, 

with 

x (0) = 1, and x'o (0) = 0. 

The solution xo(t) is just x(t, 0) given by (9.3). 

Setting to zero the coefficient of e on the leftside of (9.5), we get 

x*i +xi = -e~'cost. (9.8) 

Similarly, we get from (9.6) and (9.7) that 

*i(0) = 0. (9.9) 

xi(0) = 0. (9.10) 
The solution of (9.8) satisfying (9.9) and (9.10) is 



xi(t) = 4-cos?- -2-sinf + e~' 2sriU g 



cos? 



5 5 5 

Since both x\(t) and xo(0 are 0(1) at all times and since e is very small, we have 

kxi(f)l< \xo(t)\ 

at all times, except when t is near the points where xo(t) happens to vanish. 

One may further obtain the n th -order correction x„(t) from (9.5)-(9.7). We shall find that x„(f) is 
0(1) at all times. Thus 
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\e n x n {t)\<L le^jc-iCOI. 

The inequalities above indicate that keeping more and more terms in the series of (9.4) makes the 
approximate solution better and better. 

The series (9.4) is a perturbation series, and the method given above in obtaining a perturbation 
series is called that of regular perturbation. 

B Boundary Layer Theory 

It may be surprising, but not all problems with a small parameter can be solved by regular 
perturbation. Let us start the discussion with a simple example. 
Consider the solution of the differential equation 

f f+, = (9.11) 

with the initial condition 

y(0) = 1. 

If we set e in the differential equation above to zero, we get 

y(x) = 0. 

But this approximate solution does not satisfy the initial condition. 

In order to understand why this is so, we solve this differential equation in a closed form. We 
find that 

y(x) = e~ xl€ . 

This solution of the differential equation satisfies the initial condition for all values of e not equal to 
zero. But if we set e to zero, the solution above vanishes, and the initial condition is not satisfied. 

Indeed, considered as a function of e, the solution of the differential equation has an essential 
singularity at e = 0. It is therefore no wonder that it does not have a Maclaurin series expansion. 

As a comparison, consider the solution of the equation 

£+V-0 (9.12) 

with the same initial condition. The exact solution is 

y{x) = e~ ex . 

The approximate solution obtained by setting e in the differential equation to zero is 

y(x) = 1 

which is a good approximation of the exact solution when ex is very small. 

What makes the method of regular perturbation applicable for one but not for the other? The 
answer lies in the fact that if we set e to zero, eq. (9.12) remains a first-order differential equation, 
while eq. (9.1 1) turns into an algebraic equation. Indeed, while it is always true that the magnitude of 
ey is small compared to that of y, it is not always true that the magnitude of ey' is small compared to 
that of y. This is because, if y(x) is a rapidly varying function of x, the magnitude of the derivative of 
y(x) is much larger than y(x). In such a case, ey'(x) may be of the same order of magnitude as y(x). 
Indeed, we may easily verify that this is true for the function y(x) = e~ xk ,the solution of (9. 11). 

We encounter similar difficulties with boundary layer problems in fluid mechanics. Let us 
consider the second-order differential equation 

ey"(x) + a(x)y'(x) + b(x)y(x) = 0, (9.13) 

which generally cannot be solved in a closed form. We shall choose the interval in which (9.13) 
holds to be [0, 1] . Typically, we require that the solution satisfies the boundary conditions: 

y(0) = A,y(l) = B. (9.14) 

The parameter e is very small, while A and B are 0(1). Note that the small parameter e is the 
coefficient of y" . 

This problem is a simplified model of the boundary layer problem in fluid mechanics. As we 
know, the highest-order term in the Navier-Stokes equation in fluid mechanics is equal to vV 2 v, 
where v is the velocity of the fluid and v is the kinematic viscosity. The kinematic viscosity is often a 
small parameter, and the order of the Navier-Stokes equation decreases if we ignore the viscosity. 

Equation (9.13) can be solved with good accuracy with the use of the WKB method. (See 
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homework problem 1.) However, we gain more insights to this problem by starting afresh. 
Let us set e = 0. Then (9.13) becomes 

a(x)y' + b(x)y = 0. 

This equation is readily solved and we get 



y s (x) = exp^- j" ^j^-dx 



(9.15) 



i(x) 

The fact is that eq. (9.13) is of the second order, and has two independent solutions. But the 
approximate equation we obtain is of the first order, and has only one independent solution. Thus the 
approximate method we use yields only one of the independent solutions. This is to say that we lose 
one of the independent solutions as we set e in the differential equation to zero. It is impossible to 
satisfy the two boundary conditions of (9.14) with only one solution. 

How do we find the missing solution approximately? To get some insight into the answer, let us 
go to the much simpler case in which a and b are constants. Then the general solution of (9.13) is 
equal to 

c\e ryX + C2e r2X , 

where c\ and a are arbitrary constants and r\ and ri are the roots of the quadratic equation 

er 2 + ar + b = 0. 

This quadratic equation can be solved easily. But if we set e in this quadratic equation to zero, we get 
only one root 

r « -bla, 

which gives the solution 

£— bxia 

The solution is equal to the one denoted by y s {x) in (9.15) with a(x) and b(x) replaced by the 
constants a and b, respectively. The other root of the quadratic equation is obtained by solving the 
equation exactly. When e is very small, this root is approximately given by 

- ale. 

This root tells us that the second solution of the differential equation is approximately 

We see that this solution is a rapidly varying function of x. Indeed, it is a rapidly decreasing function 
of x if a is positive, and a rapidly increasing function of x if a is negative. 

That the missing solution is a rapidly varying function of x also explains why it is not always 
legitimate to drop the y" term in (9.13). Let us denote the rapidly varying solution by y r (x). Then 
y' r (x) is of the order of e" 1 times y r (x), and y"(x) is of the order of e" 2 times y r (x). Therefore, the 
term ey"(x) is bigger than y r (x) by a large multiple of e~ l . Again, while it is true that ey, say, is 
always much smaller than y, it is not always true that ey" is much smaller than y. 

On the other hand, y" and y' s are both of the same order of magnitude as y s . Thus it is justified to 
treat ey" as a small perturbation term. Thus we write (9.13) as 

a(x)y' s (x) + b(x)y s (x) = -ey'^x), (9.16) 

and seek the solution y s (x) in the perturbation series 

y s (x) = y?\x) + ey?\x)+~: 

Substituting this series into (9.16), we easily find that ys°\x) is given by (9.15), and that 

a{x)-^yf\x) + b{x)yf\x) = -JL.yt x \x),n = 1,2 • • .. (9.17) 

Thus ys"\x), n = 1,2, • • •, satisfies a first-order ODE with the inhomogeneous term equal to the 
negative of the second-order derivative of y.s" _1) (x). Once yi is found, we may then solve (9.17) 
with n = 1 and obtain yi l \ Similarly, all yf\x),n = 2, 3, • • •, can be obtained by successive 
iteration. 

We remark that, by treating ey" as a small perturbation, we have implicitly assumed that the 
solution is slowly varying. Since all terms in the perturbation series so obtained are slowly varying, 
the assumption is justified and the method of regular perturbation produces succesive 
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approximations of the slowly varying solution. On the other hand, this method cannot produce the 
rapidly varying solution. Clearly, if the solution is rapidly varying, it is not justified to treat ey" as a 
small perturbation term. 

To find successive approximations of the rapidly varying solution, we note that, when a(x) and 
b(x) are constants, the rapidly varying function is exp(-ax/e). Thus we surmise that, when a(x) and 
b(x) depend on x, the rapidly varying solution is approximately exp[-e _1 Ja(x)<ix]. Therefore, we 
shall put 

y r {x) = exp^-6" 1 Ja(x)Jx]v(x). (9.18) 

With the rapidly varying behavior of y r (x) taken up by the exponential factor, we expect v(x) to be 
slowly varying. Substituting (9.18) into (9.13), we get 

av' + (a' - b)v = ev". (9.19) 

Treating the term ev" as a small perturbation we may solve (9.19) with regular perturbation. We put 

v(jc) = v (x) + evi(x) +• • •, (9.20) 

and get 



as well as 



vo(*) = -4rex P rf ^rdx] (9.21) 
a(x) |_ J a \ x ) J 



a-fv n + (a' - b)v n = -frvn-u (9.22) 
ax dx 



Note that, just like (9.17) for yf\ (9.22) is a first-order differential equation for v„ and can be solved 
in a closed form. Thus we may obtain all v„ by successive iteration. We note that this prodecure 
yields us only the solution which is in the form of (9.18), where v is slowly varying. Hence it does 
not yield y s . 

If y r (x) is rapidly decreasing, we expect that it is appreciable only in a small neighborhood near 
the lower endpoint x = 0, and drops to very small values as x gets away from the origin. The graph 
for the solution (9.23) has a bump near x = 0. The solution is said to have a boundary layer near the 
lower endpoint. 

Similarly, if y r (x) is rapidly increasing, we expect it to be appreciable only near the upper 
endpoint x = 1. 

In a specific calculation, we may, instead of pluggin in a and b into the formulae given above, 
obtain both y ou t and y,-„(jc) directly from (9.13). First of all, the differential equation for y out is 
obtained by neglecting the term ey" in (9.13). The resulting equation is of the first-order and can be 
solved in a closed form. The multiplicative constant in this solution can be fixed with the boundary 
condition at the endpoint where y r (x) is negligible. The equation for y,„ can be obtained by treating 
the functions a and b as constants. 



Problem for the Reader: 
Solve approximately the equation 

ey" + (1 +x)y'+y = 0, < x < 1, 

with the boundary conditions 

y(0)=y(l) = l. 

Answer 
We have 

a(x) = 1 + x > 0, 

thus y r {x) is rapidly decreasing and there is a boundary layer near x = with a width of order e. 
For x outside of the boundary layer, we have 

i}+x)y' out + y out = 0. 

Thus we get 
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y out (x) = T ^ 7 . 

The constant c is determined from the boundary condition at x = 1 and we have 

youAx) = T | 7 . 

In particular, 

y OHf (0) = 2. 

To obtain v,„(x) near x = 0, we make the approximation 

a(x) » a(0) = 1. 

Thus 

Therefore we have 
The boundary condition 
requires that 
Thus 

We also find that 



y in (x) = 2 + c'e~ xk . 
y(0) = 1 
c' = -1. 
y ; „(x) = 2-<r* e . 



V T = ~ P~ Xl€ 

y uniform i c 

The solution is plotted in the figure below. We see that it is given by 2/(1 + x) almost everywhere. 
But the function 2/(1 + x) does not satisfy the boundary condition at x = 0. Thus the rapidly varying 
solution y r kicks in near the point x = and the solution y goes through a rapid transition inside the 
boundary layer. 

2 1 1 1 i i 



xact solution 
niform solution 



e= 0.01 




0.8 



0.2 



0.4 0.6 

x 



0.8 



Figure 9.1. 



C Turning Points 



In the preceding section, we assume that a(x) is of the same sign throughout the interval in which 
the differential equation (9.13) holds. The solutions y r {x) and y s {x), given by (9.24) and (9.15) 
respectively, are good approximations of the two independent solutions of the differential equation. 

But if xo is a simple zero of a(x), the function a(x) changes sign at xo, and the approximations 
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we have made in the preceding section fail in the neighborhood of xo- Borrowing the terminology in 
quantum mechanics, we shall call xo a turning point of eq. (9.13). 

To obtain approximate solutions near a turning point, we approximate a(x) and b(x) in (9.13) by 
algebraic functions of x, and solve the resulting equation exactly. This is the same approach as the 
one used to handle the turning point of the wave equation discussed in Chapter 7. The only 
difference is that the approximate solutions for (9.13) near a turning point are parabolic cylinder 
functions rather than Airy functions, as we shall show. 

It will be convenient get rid of the y' term in (9.13). Let us write (9.13) as 

(D + f)Dy + \y = 0. 
This tells us that, to eliminate the y' term, we put 

y = exp^-j" dx'^l^Y(x). (9.42) 

This is because, as we move the exponential factor to the left of D, we make the replacement 

Thus eq. (9.13) becomes 



D ->• D- 

2e 



(d + ^)(d-^)y + ±y = o, 



or 



2 [ 2b(x) - a (x) _ a l {x) 
I 26 4e 2 



7=0. (9.43) 



Since l/(4e 2 ) is very large, we see that an alternative way to solve eq. (9.43) is to use the WKB 
method, approximating r\{x) in (7.16) by 

\ a 1 + 2e(a' - 2b) _ \ a \ + a ' -2b 



\ 4e 2 2e 2 | a \ ' 

The WKB approximation is valid provided that a(x) does not vanish. This approximation yields 
for Y two WKB solutions, one of them a rapidly increasing function of x, while the other a rapidly 
decreasing function of x. To obtain the solutions of (9.13), we multiply these two WKB solutions by 
the exponential factor of (9.42), getting the slowly varying solution y s and the rapidly varying 
solution y r . (See homework problem 1). 

The WKB approximation fails at a turning point. Let xo be a turning point. We shall assume that 
xo is a simple zero of a(x). Then we have, near x = xo, 

c(x) « a{x - xo) 

and 

b(x) * p. 

By (7.18), the WKB approximation holds only if 

| X-X() |» V?. 

Fortunately, when x is close to xo, or more precisely 
| x - X() | « 1 , 
eq. (9.43) is approximately 



d 2 Y + [ iP-a _ a 2 (x-x ) 



2 ~1 



Y » 0, (9.44) 



dx 2 L 2e 4e 
which we will be able to solve in a closed form. 

Before we show how to do this, we first note that, when x is near xo, we have 
Thus, for x close to xo, the solution y is related to Y by 

y ~ e -a(x-x ) 2 /(4e)Y_ (9.45) 

We here take note that, if a is positive, the factor multiplying Y in (9.45) drops exponentially as 
Ix - xol increases; on the other hand, if a is negative, this factor grows exponentially as Ix - xol 
increases. 
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We shall scale the variable (x - xo) by 



v .v., = (9.46) 



Then (9.44) becomes 



where 



d 2 Y , f _P_ _ sign(a) _ 
dX 2 l\a\ 2 4 



Y = 0, (9.47) 



sign(a) = a/\a\. 

We note that, with the independent variable scaled this way, the resulting equation (9.47) is 
independent of the small parameter e. As a result, the solution of (9.47) varies by an amount 0(1) 
when X varies by an amount 0(1). Thus we say that the scale of X is unity. By (9.46), the scale of 

(x - Xo) is Je . ( We have made the implicit assumption that a is of the order of unity.) Therefore, 
there is a boundary layer with width Je near the turning point xo- 

There is a way to see directly from (9.13) that the width of the boundary layer near x = xo is of 
the order of Je . Let us choose xo = for convenience and write eq. (9.13) near the turning point as 

ey" + axy' + /fy « 0. 
To determine the scale of x, we make the change of variable 

x = SX. 



Then we have 



We see that if we choose 



d 2 _ y + aX J_ y + p y 



dX 2J dX- 



5= V?, 

the differential equation above is independent of the small parameter e. Thus both solutions 
expressed as functions of X are independent of e. This means that the scale of X for these solutions is 
unity, or the scale of x for these solutions is Je as concluded earlier. 

Equation (9.47) is a parabolic cylinder equation, which is usually written in the generic form 

+ (v + 1/2 - z 2 /4)0 = 0. (9.48) 

dz 

One of the solutions of eq. (9.48) is the parabolic cylinder function denoted by D v (z). Since equation 
(9.48) is unchanged with the replacement of z ->■ -z, D v (-z) also satisfies equation (9.48). In 
addition, (9.48) is unchanged if we make the replacements z ->■ iz as well as v ->■ -v - 1. Thus 
D- v -i (iz) and D- v -\ (-iz) are also solutions of eq. (9.48). 

We shall choose the two independent solutions of y near the turning point to be 

yi = e- ax2 '^D v ( \^x\ (9.49) 

and 

y 2 = e-"* 2/ ^D v [ - j^x ), (9.50) 





where 

_P_ _ sign(a) + 1 
led 2 

The asymptotic forms of the parabolic cylinder function are given by 

D V (X) « X v e- x2 '\ X -+ oo, 



2tt ,,„_ v _j z 2 /4 



\X\- v ~ l e x l \ X -+ -oo. 



r(-v) 

We shall replace the parabolic cylinder functions in (9.49) and (9.50) by these asymptotic forms 
when LY1 is large. 
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The qualitative behaviors of y 1 and y2 depend on the sign of a. This is because the exponential 
factor in (9.49) and (9.50) are either growing exponentially or vanishing exponentially as 1X1 
becomes large, depending on the sign of a. Thus there are two separate cases to discuss. 
Case of a > 

For a > 0, we have 

v a 

Also, the exponential function vanishes as LX1 increases. Thus we have 

(xy i+ p /a e- x2/2 , X 

J2n 



yi 



OO, 



r(l-j8/a) 



and 



y2(x) 



p/a) 



(1X1)-^°, X -> -oo 



X-P /a , X -»• oo 



(9.51) 



r(i 

-* | X | - l+ P la e~ x2/2 , X -> -oo. 
The behaviors of these two functions are schematically illustrated below: 



(9.52) 





yi 




y 2 








' — 










a<0 


a>0 


a<0 


a>0 



Figure 9.2. 

In the region x > 0, yi is the rapidly varying function, and y2 is the slowly varying solution, 
while in the region x < 0, yi is the slowly varying solution and is the rapidly varying solution. The 
role of y i interchanges with that of y 2 as one crosses the turning point. 
Case of a < 

We have 

v a ■ 

Also, the exponential function blows up as | X \ becomes large. Thus 

yi -+ X~P' a , X -»• 00, 



r(j8/a) 



| Z |-l+/?/a e X : 



/2 



X ^ -00, 



(9.53) 



and 



J2 



I2k 



-X- u P /a e x2/ \ X - 00, 



r(j8/a) 

IZI^ /a , X -00. 

We illustrate in the figure below the behaviors of yi and y2 for a < 0. 



(9.54) 
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yi 



Figure 9.3. 

In the region x > 0, yi is the rapidly varying solution and y i the slowly varying solution. And in 
the region x < 0, y\ is the rapidly varying solution and yi is the slowly varying solution. 

In the next three sections, we shall show how to match these approximate solutions valid inside 
the neighborhood of the turning point with the approximate solutions valid outside of the 
neighborhood of the turning point. 



D Turning Point at an Endpoint 



Problem for the Reader: 
Solve approximately 

ey" + 2xy' + (1 + 3x 3 )y = 0, < x < \,e <L 1, 
with the boundary conditions 

y(0) = 0, y(l) = 1. 

Answer 

Since a(x) = 2x > 0, the rapidly varying solution is a decreasing function of x. There is a 
boundary layer at x = 0. Since x = is also a turning point, the width of the boundary layer near 
x = is Je . 

The solution y r is negligible outside of the boundary layer. Thus, when x » Je , the solution is 
equal to y„ ut which satisfies 

2xy^ r + (l +3x 3 )y out = 0. 

This equation yields 

y out (x) = ^-e^'\ 

where we have made use of the boundary condition at x = 1 . 
Inside the boundary layer near x = 0, we have 

y in (x) = 

These solutions are good approximations if x « 1. 

There are two arbitrary constants, c\ and a. We determine one of these constants by matching 
Vm(x) with y„ut(x) in the overlapping region where both approximations are valid. The parabolic 
cylinder function solution is valid when x « 1, while by (7.18), the WKB solutions are valid when 




d e 



\« 1, 



dx x 
or 

x » y?. 

Thus the overlapping region is 

y? « x « 1. 

When x » y? we have 



y in (x) *c 2 (26) 1/4 x 



1/4^-1/2 
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Also, when x « 1, we have 

youtix) « x~ m . 

We note that both v,„(x) and y ou t(x) are equal to a constant times x~ m in the region 1 » x » Je . 
This is an indication that this region is indeed the overlapping region in which both approximations 
hold. Joining y ,•„(*) with y out in this region, we get 

c 2 = (2e)~ m . 

From the boundary condition at x = 0, we get 

c\ = -C2- 

Thus we have 



1/4 



-D. 



1/2 



f -V ) +£>-l/2l - 



2. 



As a final observation, we note that y OHt (0) is infinite. But as we continue y ou t{x) into the region of 
the boundary layer, it turns into 



r7(2e) 



(2ey" 4 D- m (-JJxy 



At x = 0, the expression above is equal to (2e) 1/4 Z)_i/2(0), which is a large number but is not 
infinity. 
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